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New boundary conditions are constructed and tested numerically for a general first-order form 
of the Einstein evolution system. These conditions prevent constraint violations from entering the 
computational domain through timelike boundaries, allow the simulation of isolated systems by 
preventing physical gravitational waves from entering the computational domain, and are designed 
to be compatible with the fixed-gauge evolutions used here. These new boundary conditions are 
shown to be effective in limiting the growth of constraints in 3D non-linear numerical evolutions of 
dynamical black-hole spacetimes. 



I. INTRODUCTION 

The Einstein system can be written as a set of evolu- 
tion equations that determine how the dynamical fields 
change with time, plus constraint equations that must be 
satisfied by the physically relevant field configurations. 
The evolution equations ensure that the constraints will 
be satisfied within the domain of dependence of the ini- 
tial data if they are satisfied initially. But this does not 
guarantee that constraints that are initially small (rather 
than precisely zero) will remain small, or that constraint 
violations will not enter a domain through its timelike 
boundaries. Indeed, the rapid growth of constraint vi- 
olations from small (truncation or even roundoff level) 
values in the initial data continues to be one of the ma- 
jor problems for the numerical relativity community. 

Constraint violations in the continuum evolution equa- 
tions (as distinct from their discrete numerical approxi- 
mations) have at least two different causes: a) bulk con- 
straint violations, in which existing violations are ampli- 
fied by the evolution equations, and h) boundary viola- 
tions, in which constraint violations flow into the domain 
through timelike boundaries. A variety of techniques 
have been introduced to control bulk constraint viola- 
tions in numerical solutions of constrained evolution sys- 
tems; these include constraint projection (where 
the constraint equations are re-solved whenever the con- 
straints become too large), fully constrained evolution 3, 
1 S i, 1 E im (where some dynamical fields are 
determined at each time step by solving the constraint 
equations rather than their evolution equations), and dy- 
namical constraint control [T^ IT^ IT^ (in which the evo- 
lution equations are modified dynamically in a way that 
minimizes the growth of the constraints). Techniques 
have also been proposed to control boundary constraint 
violations; these include the construction of special forms 
of the evolution equations that prevent constraint viola- 
tions from flowing into the domain |l5l Il6| , and special 
boundary conditions that prevent the influx of constraint 
violations in more general forms of the evolution equa- 
tions II 111 111 111 111 S 111 111 111 111 111 111 



Recent studies have shown that bulk constraint viola- 
tions are not effectively controlled by constraint projec- 
tion Pi or dynamical constraint control 'l2'| techniques, 
unless the influx of boundary constraint violations is also 
controlled separately. Thus the explicit control of bound- 
ary constraint violations appears to be an essential re- 
quirement for effective constraint control in the Einstein 
system. The primary purpose of this paper is to con- 
struct and then test numerically the special boundary 
conditions needed to prevent the influx of constraint vi- 
olations in one (rather general) first-order form of the 
Einstein evolution system: the multi- par ameter general- 
ization of the Frittelli-Reula system l2a introduced by 
Kidder, Scheel, and Teukolsky (KST) |29l|. These con- 
straint preserving boundary conditions are constructed 
here using the strategy first outlined by Stewart jl^ for 
the Einstein equations. The idea is to decompose the 
constraint fields into incoming and outgoing parts, based 
on the characteristic decomposition of the constraint evo- 
lution equations. The incoming constraint fields are con- 
trolled {e.g. set to zero) at each boundary point, and 
these conditions then serve as boundary conditions for 
the principal dynamical fields of the system. The con- 
straint preserving boundary conditions derived here gen- 
eralize earlier work by being applicable to generic non- 
linear field co nfig urations 17, 20], having no symmetry 
requirements |l8| . and allowing arbitrary values of the 
gauge fields {i.e., the lapse and shift) 0, |12- These 
new boundary conditions are also tested here under more 
challenging conditions — non-linear 3D evolutions of dy- 
namical black-hole spacetimes — than previously consid- 
ered. 

A secondary purpose of this paper is to introduce and 
test boundary conditions for the physical (gravitational 
wave) degrees of freedom of the Einstein system. These 
new physical boundary conditions control the influx of 
the radiative part of the Weyl tensor, and generalize (to 
the generic 3D case) the boundary conditions of this type 
proposed by Bardeen and Buchman [T^ . Boundary con- 
ditions are also introduced and tested here for those dy- 
namical fields that are not fixed by the physical or the 
constraint preserving boundary conditions. The bound- 
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ary conditions for these extra "gauge" degrees of freedom 
are set in a way that is consistent with the "fixed-gauge" 
{i.e., time independent lapse density and shift) evolutions 
used in the numerical tests described here. 

The remainder of this paper is organized as follows. 
The basic KST form of the Einstein evolution sys- 
tem is reviewed in Sec. ^ and additional technical de- 
tails (which considerably generalize previously published 
work) on the characteristic decomposition of the dynam- 
ical fields of this system are given in Appendix El The 
constraint evolution equations associated with the KST 
system are described in Sec. IIIII and additional details 
on the hyperbolicity of the constraint system are given in 
Appendix The principal new analytical results of this 
paper arc contained in Sec. lIVI where constraint preserv- 
ing boundary conditions are derived for the KST system. 
Two types of constraint preserving boundary conditions 
are presented here: the first, in Sec. IIV Al are based on 
the characteristic constraint decomposition ideas of Stew- 
p^ . while the second, 



art [19|J, while the second, in Sec. IIV 51 are generaliza- 
tions of a simpler (but less generally applicable) type of 
constraint preserving boundary condition which has been 
used effectively for the scalar-wave system Q. Physical 
boundary conditions that control the influx of the ra- 
diative part of the Weyl tensor are presented in Sec. ^ 
Boundary conditions for the remaining "gauge" dynami- 
cal fields are given in Scc. lVIl We test the efficacy of these 
new boundary conditions by performing numerical evo- 
lutions of various perturbed and unperturbed black-hole 
spacetimes. The basic numerical methods used in these 
tests are described in Sec. I VIII The tests themselves are 
described in Sec. IVIIII tests performed on unperturbed 
black holes are given in Sec. IVIII Al tests on perturbed 
black holes are presented in Sec. lVIIIBl and finally a mild 
"angular" numerical instability that appears in some of 
these tests is discussed in Sec. lVIirO The major results 
of this paper are summarized and outstanding questions 
raised by this work are discussed in Sec. IIXI 



II. PRINCIPAL EVOLUTION SYSTEM 

The form of the Einstein evolution system used here 
is a first-order formulation introduced by Kidder, Scheel, 
and Teukolsky (KST) , which generalizes several ear- 
lier forms of the Einstein system '2s , 'so', 'si'l . This system 
consists of first-order evolution equations for the spatial 
metric gij, the extrinsic curvature Kij, and the spatial 
derivatives of the metric D^ij = d^gij/'i. The principal 
(or highest derivative) parts of these evolution equations 
are given by 



dtD 



^tg^J ~ N"^ng^J, (1) 

\l + 27o)g^'*<5"(,(5^,) 
-(1 + 72)5"'*5^.<5=,) - (1 - l2)g''S\,5'',^ 
+g'''S-,S''j + 27iff"['5""5z,l dnD,,^, (2) 
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dnKbc, 



(3) 



where ~ indicates that terms algebraic in the fields are 
not shown explicitly. This form of the equations assumes 
that the lapse density. 



Q = log {Ng- 



(4) 



and the shift N'' are specified a priori as functions of 
the coordinates rather than being evolved as indepen- 
dent dynamical fields. The parameter^ 70 that appears 
in these equations is part of the definition of the lapse 
density, Eq. (0J, while the parameters 71, 72, 73, and 74 
were introduced by adding multiples of the constraints 
to the standard ADM form of the evolution equations 
(see KST |2^). This form of the Einstein equations is a 
quasi-linear first-order system that can be written more 
abstractly as 



dtu" 



pa 



(5) 



where u" = {gij , Kij , Dj^ij } is the thirty dimensional vec- 
tor of dynamical fields, and the quantities A^^^fj and 
depend on but not its derivatives dku". In this paper 
the Greek indices a and (3 are used to label elements of 
this space of dynamical fields. 

Boundary conditions for hyperbolic evolution systems 
like Eq. ((SJ are imposed on the incoming characteristic 
fields of the system at each boundary point. These char- 
acteristic fields are defined as follows. Given a direction 
field rik {e.g. the outward directed unit vector normal to 
the boundary) we define the left eigenvectors e"a of the 
characteristic matrix UkA'^'^p by 



(6) 



where is the eigenvalue (also called the characteristic 
speed). Greek indices with hats, like d, label the various 
eigenvectors and corresponding eigenvalues of the sys- 
tem. There is no summation over d in Eq. ©. These 
eigenvectors form a complete set in any strongly hyper- 
bolic system of equations, so the matrix e"a is invertible 
in this case. The projections of the dynamical fields u" 
onto these characteristic eigenvectors are called the char- 
acteristic fields of the system: 



(7) 



Boundary conditions must be imposed on any character- 
istic field having a negative characteristic speed < 
{i.e., an incoming field) at a particular boundary point. 



^ The pa ram eters 70 to 74 used here are related to the parameters 
of Ref. |2E| by 70 = it, 71 = 7, 72 = 73 = rj, and 74 = x- 
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The characteristic fields of the KST system are the 
cohection of fields u" = {9^J,Z\ Z}, Zf, Z*, Z^, Zf^^-, 

f/f=^, C/3±, C///}. The definitions and exphcit ex- 
pressions for these fields are given in Appendix 1X1 The 
characteristic fields C/^*, and U^'^ have characteris- 
tic speeds iwi , ±i;2 , and ±z)3 respectively relative to the 
hypersurface orthogonal observers, where 



III. CONSTRAINT EVOLUTION SYSTEM 

Boundary conditions capable of preventing the influx 
of constraint violations cannot be formulated without a 
complete understanding of how the constraints propa- 
gate. The constraints associated with the (vacuum) Ein- 
stein evolution system are 



270, 

-73(1 - 372 - 470) 



-74(1 + 670), 



1/ N 1 

-(l-f 27i)(2 - 73 274) - -7273- 



(8) 
(9) 

(10) 



The characteristic fields XJ^p have speeds ±1, and {5^^, 
, Zl, Zf, Zf, Zfj, Zf-j} have characteristic speed zero 
relative to the hypersurface orthogonal observers. 



The KST evolution system is strongly hyperbolic if and 
only if the characteristic fields are linearly independent. 
This is the case if vi ^ 0, V2 ^ 0, V3 ^ 0, and vi ^ V3. 
The system is also strongly hyperbolic for vi — V3 ^ 
and 1 -I- 3vf — 4w|. In the strongly hyperbolic case the 
fundamental characteristic fields have thirty independent 
components: the spatial metric gij (with six independent 
components), five scalars {Z^, U^^, U^^}, five transverse 
(to n'') vectors {Zf, Zf,Zf, Uf^} (which have a total of 
ten independent components), three transverse-traceless 
second-rank tensors {Z^^, U^^} (which have a total of 

is symmet- 

and one transverse-traceless third- 



1-0 

seven independent components, since Uf^ 
ric while Zf, is not) 



rank tensor Z^^^- (which has a total of two independent 
components). The KST system also admits a positive- 
definite symmetrizer matrix Sap (and so is symmetric 
hyperbolic) under fairly weak conditions on the KST pa- 
rameters, which are discussed in detail in Appendix B of 
Ref. [33. 

The characteristic speeds relative to the coordinate 
frame, i.e., the eigenvalues of Eq. ©, are ±Nvi — 
UkN^, ±Nv2 ~ UkN^, ±Nv3 - UkN^ and ±N - nuN^ 
for the fields C/i±, t/^^, and J/// respectively. And 
the speeds are —nkN^ for the "zero-speed" fields {gij, 
Z\ Zf, Zf, Zf, Z^., Zf^^.}. Boundary conditions must 
be supplied for each incoming characteristic field, i.e., for 
each field whose coordinate characteristic speed is nega- 
tive at a particular boundary point. The fields , , 
, and Ufj~ will be ingoing at most timelike bound- 
aries; the zero-speed fields {g.^, Z^, Zf, Zf, Zf, Zf^, 
^kij } ^ill be ingoing at any boundary where the shift N'' 
is directed out of the computational domain, i.e., when- 
evcTUkN'' > 0. Boundary conditions must be formulated 
therefore for each of the fields {gij, Z^ 
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c 

Ckij 



klij 



dkgij — 2Dkij, 



(11) 

(12) 

(13) 
(14) 



Here C is the Hamiltonian constraint, Ci is the momentum 
constraint, while Ckij f^nd Cknj are auxiliary constraints 
associated with the introduction of the dynamical field 
Dkij needed to make the KST system first order. All of 
these constraint fields are zero for the physical solutions 
of the (vacuum) Einstein evolution system. Note that 
Ckij and Ckiij are not completely independent: there is 
a second-class constraint Ckiij = 9[;Cfc]y. However, Ckij 
and Ckiij must both be retained and treated as indepen- 
dent in order to write the evolution of the constraints as 
a closed first-order hyperbolic system, which is the goal 
of this section. 

The evolution of the constraints is completely deter- 
mined by the evolution of the dynamical fields of the 
principal evolution system. Using the KST evolution 
equations, Eqs. Q JSJ, and the definitions of the con- 
straints, Eqs. Hll|) (|14|l . it is straightforward to show that 
the principal parts of the constraint evolution equations 
are 



dtC ~ N''dkC-^{2-^3 + 2^^)Ng'^d,C,, 

dtC, ~ N''dkC,~il + 2-fi)Nd,C 
1 



(15) 



+ -Ng''g-'' (1 - j2)dkClaM + (1 + l2)^kCa^lt 



-(1 -I- 2jo)dkCliab 



dfCkij 
dtCkiij 



N'^dkC. 



N-^d^CkHj + ^l3N{g^[idk]C, + g^[idk]C,) 
+l4Ngijd[kCi]. 



(16) 
(17) 

I 

(18) 



More abstractly, the constraints satisfy a quasi-linear 
evolution system of the form 



dtc^ 



A'^^Bdkc'' 



t bc 



(19) 



where — {C, Ci, Ckij, CkUj} are the constraints defined 
in Eqs. (|11|I - H14|I . A^"^ b depends on the dynamical fields 
m", and B depends on and dku". We use upper 
case Latin indices such as A and b to label the constraint 
fields. If the constraint evolution system is hyperbolic, 
then Eq. H19|l guarantees that the constraints will vanish 
everywhere if they are zero at some initial time, and if 
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boundary conditions are chosen to prevent the influx of 
constraints through timehke boundaries. 

In order to determine whether this constraint evolu- 
tion system is hyperbohc, we evaluate the characteristic 
fields associated with Eq. ((T^ . It is straightforward to 
show that the fields = {Ck^J, Z\ Zf, Zf, Zl°, Z^l;, 
Uf^, U^^} are characteristic constraint fields, where the 
individual components of are defined by 

Z' = 73C-(2 - 73 + 274Kn'C^;, (20) 

Zf = n''P',[j4Cl-{j3 + ij4)n^n''Ckiat], (21) 

Zf = n^P\{iCli + 2Cli-7n'^n'-Ckiab), (22) 

4" = P\P'flu (23) 

Zll = {P\P',^\P.,P''')CI„ (24) 

(25) 

C/f ± = ±n'=P', [(1 + 27o)Ci; + 2C|;] - 272C(2,;)] 

+2v2P\Cu (26) 
f/6± = (l + 27i)C±f3n'=Cfc-72nVC|,. (27) 

The V2 and 173 that appear in Eqs. H20() - (|27|l are given by 
Eqs. ® and and C}^ and C^- are defined by 

Cjj = g^^Cijku (28) 
Cl = g'''Ck^,l. (29) 



IV. CONSTRAINT PRESERVING BOUNDARY 
CONDITIONS 

The constraint characteristic fields and are 
incoming fields at most timelike boundaries, and the 
fields {Cfe,„ Z\ Zf, Zl Z}f, Zfl, Zlf^i} are incom- 

ing at any boundary where n^N^ > 0. Therefore we 
must choose boundary conditions on the incoming char- 
acteristic fields of the principal evolution system, u", in 
a way that controls these incoming constraint charac- 
teristic fields, c^. We use two different approaches to 
accomplish this for the KST system. The first approach 
(introduced by Stewart [I^ and developed by Calabrese, 
et al. re-expresses the incoming characteristic con- 
straint fields, , in terms of the principal characteristic 
fields, u". This results in a set of Neumann-like bound- 
ary conditions on certain incoming characteristic fields 
u". These Type I boundary conditions are discussed 
in more detail in Sec. IIV Al The second (less general) 
approach (introduced by Hoist, et al. 0) uses a more 
direct Dirichlet-like boundary condition for certain in- 
coming characteristic fields . These Type II boundary 
conditions are discussed in more detail in Sec. IIV Bl 

We note that it is not possible to use these methods 
to derive boundary condition for all of the characteristic 
fields (of the principal system) that need them. In par- 
ticular it is not possible to obtain boundary conditions 
for the fields Zf, , and Ufi~ in this way. The bound- 
ary conditions for these fields are determined by physical 
and gauge considerations, which are discussed in detail 
in Sees. El and ED 



The characteristic constraint fields J7f and U have 
characteristic speeds ±V2 and ±^3 respectively, while 
the characteristic constraint fields {Ckij, Z"^ , Zf, Zf , 
Z}^ , Zfj, Zj'^f^i} have characteristic speed zero relative 
to the hypersurface orthogonal observers. The char- 
acteristic constraint fields have forty independent com- 
ponents: Ckij (which has eighteen independent compo- 
nents), three scalars {Z'^ ,U^^}, four transverse (to n*^) 
vectors {Zf,Zf,U^^} (which have a total of eight in- 
dependent components), one antisymmetric transverse 
second-rank tensor Zjj' (which has one independent com- 
ponent), one transverse traceless second-rank tensor Zf^ 
(which has three independent components), and one to- 
tally trace-free fourth-rank tensor ^^j,; (which is anti- 
symmetric in its first two indices, symmetric in its last 
two indices, and so has seven independent components). 
These characteristic constraint fields are linearly inde- 
pendent so long as f 2 7^ and V3 0. Thus the charac- 
teristic constraint evolution equations are strongly hyper- 
bolic whenever the principal evolution system is strongly 
hyperbolic. In Appendix IbI we summarize the conditions 
under which the constraint evolution system is also sym- 
metric hyperbolic. 



A. Type I Boundary Conditions 

The incoming constraint characteristic fields Ckij, Z'^ , 
Zf, Zf, Z\], Z\%„ C/f-, and C/^" all depend on nor- 
mal (i.e., perpendicular to the boundary) derivatives of 
the principal characteristic fields m". A straightforward 
but lengthy calculation shows that these characteristic 
constraint fields can be expressed as follows: 

n^Cuj = d^gij - In^Dkij, (30) 
= -d^Z^ - (2 - 73 + 2-fi)P''''n^n''daD,ab 

-2j3P^'g^^daD[,,], + ^^.Da.^iSD^"^ - 2/?-"), (31) 

Zf = d^Zf-lj^g"' - (73 + 3j4)n''n'']n'P'',ddD,ab, (32) 
Zf = d^Zf - [(3g''^ - 7n-n')P^ - 2P\P"']n-ddD,^b, 

(33) 

= dj_Zfj 

-{P\p^j - iPyP'^'')(.g^'^5ai^eM - P'%Dabd), (34) 
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"'^'feS'^icafc = d^Zl, - n^P^r'ddD^at, (35) 

+2v2{g'''g'"'P''^ + g^'g^'P^ - 2g''-g''' P'',)D,abK,k 
- [(1 + 72)n".g'" - (1 + 27o)n'5™"] P\dkDi^n 
-[(1 - 72)"'P"\; - (1 + 72)"'"P'.]P'^'"5feA™„, (36) 

= - v^n'P^'^dkK.i - 72P'="n'n"afcA„m 

+ (1 + 27i)(2g"['P"l'=afcAm„ + 5'''-'5""'^.j--?fah) 
-«3(ff'^^5''n^ + g'^g'^'n'' - 2g-'^g'"n'')D,abK,k 

-i(l + 27i)(5^^5*^g«'' + 2g'='^5*''.g^'^ + Sg'^I^ff'^l^S^'' 

-ig'''g"'g^'')Dk,jD,ab, (37) 

where P^"'' is defined by Eq. lAlSI) . and the quantities 
d^gi-j, d±Z^, d±Zf, d±Zf, d±Zfj, d±Z^^^, d±_Uf^ and 
d±U^~ are components of d±u°', the characteristic pro- 
jections of the normal derivatives of the dynamical fields. 
These projected normal derivatives are defined by 

d^u" = e^pn^dku'^. (38) 

We note that Eqs. H3U|) -- (|37|I express these incoming char- 
acteristic constraint fields, c"^, in terms of the normal 
derivative of one of the principal characteristic fields, 

, plus terms that depend on derivatives tangent to the 
boundary and terms that are algebraic in the fields. 

The simplest Type I constraint preserving boundary 
conditions are obtained by setting the incoming compo- 
nents of the characteristic constraint fields to zero. Us- 
ing Eqs. (|Sn|l - l|T7|l . these expressions become boundary 
conditions for the normal derivatives of the incoming dy- 
namical fields d^gij, d±Z^, d±Zf, d±Zf, d±Zf^, d^Z^^j, 

di_Uf- and d^U^--. 

d^g.j ^2n''Duj, (39) 
d^Z^ - -(2 - 73 -f 2^i)P-''n-n''daD,db 

+ ^73 (SD^ab]'Dc''' - Dab'D''-, - KabK-' + K^) 

-2-fsP-'g"'daD[,,]^ + ^^sDabciSD-"' - 2i?^°'), (40) 

d^Zf = [74.9'^" - (73 + 374)n'^7i^] n-P^^ddD^ab, (41) 
d^Zf = [{Sg^" ~~ ln'^n'')P\ ~ 2P'^,P'"^] jfdiD,,^, (42) 

d^Z% = {P\P'j - ^P^jP^''){g'"'daDM " P'^d.Dabd), 

(43) 

d^Zl^ = n-PtfddD,,,, (44) 

d^Uf- = 2v2{P'''P\ - g^'P^)dkK^i 

+2«2(5"^5''^'. +5''5"''^'. - 2g'''g''^P\)D,abK,k 
- [(1 + 72)n"ff'" - (1 + 27o)ri'g""] P\dkDi„,n 
-[(1 - 72)n'P"% - (1 + 72)n"P',]P^-"afcA™„, (45) 



d±(73- = v^iu'P^'^dkKji + 72P'-'"ri'"'"5feA™„ 
-(1 + 27i)(2g"['P"l'=afeA™„ + g'^'g^^^^K.^Kab) 
-^«3(5'^^ff''n^ + g'^g'^^n'^ - 2g--g''^n'')D,abK,k 

+ 1(1 + 271) (/^g^^g"" + 2/-5*''.g^'^ -I- Sg'^^^g^^^ g^" 

--3g''^g'''g^'')Dk^JD,ab■ (46) 

More general Type I boundary conditions of this type 
can also be obtained by setting the incoming constraint 
characteristic fields to multiples of the corresponding out- 
going fields: i.e., by setting Uf~ = fJ-sUf'^ and U^" — 
fiQU^^. These boundary conditions can therefore be gen- 
eralized by adding the terms fJ-sUf^ and —fiaU^^ to the 
rights sides of Eqs. 145|) and (|46|l respectively. 

We point out that the combination of constraints, 

Zf - IbP^.N^'n'^n^Zlla = 16P'^[^P'^U'arfi?eah, (47) 

involves no normal derivatives of D^ab- Thus, arbitrary 
multiples of this combination of constraints could be 
added to the Type I boundary conditions for , Z^ and 
Uf^ in Eqs. gTJ, (g^l and 103) without changing the 
basic structure of these conditions at all. We find that 
the addition of these terms does change the stability of 
numerical evolutions. But at present we have not found 
a systematic way to optimize the addition of these extra 
constraint fields (short of brute force numerical testing), 
so the numerical evolutions presented here use the sim- 
plest analytical forms given in Eqs. (|41|l . (|42|l and H45|) . 

We find it convenient to impose these Type I con- 
straint preserving boundary conditions via the Bj0rhus 
method 1331 (which we commonly use in our numerical 
code) in the following way. The characteristic projection 
of the principal evolution system can be written as 

dtu" + = e%(-A'=\P'fea,u" + ^^5), (48) 

where dtw" represents the characteristic projection of the 
time derivative of the dynamical field: 

dtu" = e^'pdtuP. (49) 

The terms on the right side of Eq. (|48|l depend on the 
boundary values of the dynamical fields and the deriva- 
tives of these fields tangent to the boundary surface. We 
leave these terms unchanged. However, we replace the 
term d^u" on the left side of Eq. (|48|l by its desired 
value dx_u°^^. For the fields 5y, Z^ , Zf, Zf, Zf^-, U^' 

and C/^", dx^u^^ is given by Eqs. 

This method of imposing these boundary conditions 
can be implemented numerically in a simple and elegant 
way. Let us define the quantity 

Au" = e"0(-v4'=^„afe?/" + Pf"), (50) 

which is just the characteristic projection of the right side 
of the full evolution equation, without any modifications 
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for boundary conditions. At the boundary, Eq. (|48|) then 
takes the form 

dtw" = AM"+t'(a)(d-L""-rfj-"y- (51) 

In this expression, d±u^^ represents the desired values 

of the tangential derivatives at the boundary, and d±u°' 
denotes the actual values of the tangential derivatives as 
defined in Eq. (|38|l . When d±u°' is expressed in terms 
of the constraints on the boundary using Eqs. 130() - (|37|l . 
and when d±u^^ is determined from Eqs. (|39|l -14t) |l . then 

the term d±u°' — d^u"^^ simplifies considerably. Thus in 
particular, the boundary conditions for g^, Z^, Zf, Zf, 
Zfji ^feii i Uf" and U^' can be written as 



dtgij 


= Dtgij - 






(52) 


dtZ' 


= + 


nkN''Z\ 




(53) 


dtZf 


= A^f- 


nuN^Zl 




(54) 


dtZf 




UkN^Zl 




(55) 


dtZ% 


= DtZl - 


^nuN^Zll 




(56) 


dtZlij 


= D^Zl^ 


n lum d pcab y 12 




(57) 


dtUf- 


= AC/f- 


+ {Nv2 + nkN''){Ut 




(58) 


dtU""- 




- {Nv^ + nkN*'){U^- 




-). (59) 



The Dtu°' that appear in these expressions are to be eval- 
uated using Eq. (|50() . while the constraint fields Ckij, Z'^ , 
Zf, Zf, Zfl, Zll^, and C/^" are to be evaluated nu- 
merically using their definitions in Eqs. p3|l . H2Uf) . (|21|1 . 
(EH, (Ell), HH), lEni) and (|2ZIl respectively. 



B. Type II Boundary Conditions 

A second type of constraint preserving boundary con- 
dition can be imposed on those characteristic fields that 
represent various projections of the field P'^iDkij- The 
infiux of constraint violations that could cause P^iDkij 
to differ from P^idk.gij /"i. can be prevented by enforc- 
ing the equality of these fields on the relevant bound- 
aries. This is a Dirichlet-like boundary condition that 
sets P'^iDkij — P'^idkgij/2. This type of constraint pre- 
serving boundary condition has been used successfully in 
the scalar wave system by Hoist, et al. T|. The char- 
acteristic fields zf, zf, Zfj, and Z^^j are composed en- 
tirely of various components of P^iD^ij, which represent 
the derivatives of gij that are tangent to the boundary. 
These Type II boundary conditions can be imposed using 
the Bj0rhus method by setting 

dtZf = i[74P'^'-(73 + 274)n'^n'']F^,a,at5a6, (60) 

dtZf = i [{ZP'^' ~ ^n'^n')P^, - 2P\P'-]dAgabm 

dtZl = ^{2P^P'^, ~ P^.P-'^ydcdtgab, (62) 

dtZl, = \p'kt'dA9ab. (63) 



at any boundary where these fields are incoming. The 
quantity dcdtgab that appears on the right sides of 
Eqs. H6()|l -H63 |l is to be evaluated by taking the indi- 
cated tangential spatial derivatives of the right side of 
Eq. H52I) . These expressions are obtained using the ex- 
pressions for Zf, Zf, Zfj, and Z^-j from Eqs. l|A2p . 
(IA3(I . (|A5|I . and (|A6(I . and replacing the terms dtDcab 
by dcdtgab/2. We note that these fields Zf, Zf, Zf^ and 
can also be controlled using the Type I boundary 
conditions described in Sec. IIV Al We have tested these 
Type II boundary conditions numerically, and find the 
results to be essentially identical to the results described 
in Sec. IVIIll for tests with Type I boundary conditions. 

V. PHYSICAL BOUNDARY CONDITIONS 

In this section we derive boundary conditions for the 
physical components of the dynamical fields, Uf~ , which 
represent the real gravitational wave degrees of freedom 
of the system, at least asymptotically. Our strategy is 
based on the proposal of Bardeen and Buchman > and 
similar ideas used by Reula and Sarbach in the con- 
text of the Maxwell system. These ideas are adapted to 
the general 3D Einstein system by analyzing the dynam- 
ics of the Weyl tensor as determined by the Bianchi iden- 
tities. We then infer boundary conditions on the dynam- 
ical fields of the KST system that control the incoming 
radiative parts of the Weyl tensor. We begin by following 
the usual practice of decomposing the Weyl tensor into 
electric and magnetic parts, 

E^, = C^aurT^T^ (64) 

B^i, — — C^cjcTT- e i/pT , (65) 

where T^^ represents a timelike unit vector. [In this pa- 
per letters from the latter part of the Greek alphabet (/x, 
V, etc.) represent four-dimensional spacetime coordinate 
indices.] In our analysis here we assume that is or- 
thogonal to the t ^constant spacelike hypersurfaces of 
the standard 3-1-1 decomposition of the geometry. In this 
case we can express the electric and magnetic parts of 
the Weyl tensor in terms of standard 3-1-1 quantities: 

E,, = R,,+KK,,~K,''Kkj, (66) 

= {VkKu)e''',, (67) 

where Rij, Kij, Vi and e^fc represent the three- 
dimensional Ricci tensor, the extrinsic curvature, the 3D 
spatial covariant derivative, and the 3D totally antisym- 
metric tensor respectively. These electric and magnetic 
parts of the Weyl tensor are symmetric and traceless in 
any vacuum spacetime. 

The Weyl tensor satisfies the Bianchi identities, 

= ^laCru.]^,, (68) 

in any vacuum spacetime, where here represents the 
4D covariant derivative, and these equations imply a sys- 
tem of evolution equations for the Weyl tensor. The 3-1-1 
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representation of the principal parts of these evolution 
equations can be written ||35j| in the form 



~N[dkB, 



kl 



(69) 
(70) 



which is reminiscent of the Maxwell system. The charac- 
teristic fields of this vacuum Weyl tensor evolution sys- 
tem are 



7l3 



7± 



8± 



n''n''Eab, 

1 



(71) 
(72) 
(73) 



using 



Eq. (EU: 



X l^n'daKbc + 2j2P'''d[aD,]bd + l2n'n''daD,bd 

+ 2P'''dcD^da]b + 2g'=''daD[,,]a - KKab + KacK% 

+ {2D''', - D^',){2Dabd - Ddab) + AD'\Dydc\b 
-Da'^Dbcd + 2n''K\,D,^t,a + ^n^ D^^^^^.K^f 

(76) 

In deriving this expression, we have used the fact that 
the constraint characteristic field Zl^ can be expressed 
in terms of using Eq. Setting Uff = C/y~lbc 

(and = 0) gives us the desired Neumann-like bound- 



(P(^P^), - -P^.P'^'') {Eab T ea'^naB^b) , (74) ary condition for the principal characteristic field Uf. 



where n is a spatial unit vector. The characteristic fields 



Z and Z have characteristic speed zero, Uj^ have 
speeds ±1/2, and Ufj^ have speeds ±1 relative to the 
hypersurface orthogonal observers. These fields propa- 
gate with coordinate speeds —UkN^, zLN/2 — nkN'', and 
±7V — rikN^ respectively. An equivalent expression for 

Uft^ in terms of the standard 3-1-1 variables is 

''J 



+2^2P''-d^aDc]bd + l2n^n''daD,db + 2P--'d,D^aa]b 
+2g"'daD^bc]d - KKab + KacK\ + AD'\Dydc]b 

+ {2D-'\ - D''\){2Dabd - Ddab) - Da^'^Dbcd 



+2n''K\aDc]bd + ^n'^D^bdWcKa] 



(77) 



8± 



{P\P\, 



^Py P"'') (Pab + KKab ~ Ka'K,b 

Tn''V,Kab±n-V(aKb)^- (75) 



and C/j^ are pro- 



The vacuum characteristic fields 

portional to the Newman-Penrose [33 components of the 
Weyl tensor ^'4 and ^I^o, respectively. 



The true gravitational wave degrees of freedom are rep- 
resented by the fields U^j^. The best (gauge invariant) 
way to impose boundary conditions on the physical de- 
grees of freedom of the Einstein system is therefore to 
require that the incoming part of this field has prescribed 
values: U^^ ~ ^i^" Ibc boundaries. This condition 

is equivalent to fixing the Newman-Penrose component 
^0 = ^olbc boundary. Since U^r falls to zero as 

in an asymptotically flat spacetime j36j |. it is not un- 
reasonable simply to set Uf^ = as a physical boundary 

condition. However, since 7^ in the Kerr geometry, 
we think it is more consistent to freeze U^' to its initial 
value by setting C/fj~lbc = ^fj"(^ = 0). 

In order to write the boundary condition on Uf^ as 
a condition on the principal dynamical fields, we first 
express in terms of the principal characteristic fields 



This physical boundary condition can be imposed nu- 
merically with the same technique used for the Type I 
constraint-preserving boundary conditions discussed in 
Sec. IIV Al In particular we set 

dtUtf = DtUtf 

-{N + UkN^) [Uf- - Uf- + 724^,)] . (78) 

The term DtUf~ is to be evaluated using Eq. ISUJ), 
while ZlJ- and are to be evaluated numerically using 
Eqs. (Pl|l and ((75|l respectively. 



VI. GAUGE FIXING BOUNDARY 
CONDITIONS 



Next we turn our attention to finding boundary condi- 
tions for the two fields and Zf that were not fixed 
by either the constraint preserving boundary conditions 
in Sec. lIVI or the physical boundary conditions in Sec.lVl 
Since these fields are not fixed by the constraints or by 
physical considerations, they must be gauge fields in ef- 
fect. Boundary conditions on and Zf should be cho- 
sen in a way that is consistent with the gauge conditions. 
In the evolution system used here, we assume that the 
lapse density and the shift are frozen to their initial val- 
ues. One reasonable boundary condition for these gauge 
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fields then is simply to freeze them at their initial values 
on the boundaries: 



d,Zf 



0, 
0. 



(79) 
(80) 



We have also found that a useful boundary condition 
for can be obtained from one of the standard equa- 
tions used to fix the lapse: dtK = 0. Expressing dtK in 
terms of the time derivatives of the principal character- 
istic fields, we find 

dtK = g'^dtK,j~K'^dtg,„ 



3-q 
4w3 



3- 



+ dtU^-)-K'=^tg^J. (81) 



[The parameter q that appears in this expression is the 
ratio of characteristic speeds defined in Eq. (|A14|) ]. Set- 
ting dtK — Q therefore provides a boundary condition for 
the incoming characteristic field U^~: 



dtU'- - -d,[/l+-!^li^-^(d,C/3+ + rf,[/3 

+AvlK'^^tg^J. 



(82) 



The quantities dtU^^ and dtgij on the right side of 
Eq. (|82|l are evaluated with the appropriate boundary ex- 
pressions for these quantities. We find that using Eq. 
in evolutions of black-hole spacetimes is fairly effective 
in controlling the growth of perturbations in with 
low spherical harmonic indices. But at the same time 
using this boundary condition causes unstable (and non- 
convergent) growth of perturbations with large spheri- 
cal harmonic indices. We solve this problem numerically 
by applying this boundary condition only after applying 
a filter that removes all spherical harmonic components 
above £ — 2 from Eq. (|82|1 at the boundary. Thus we use 
Eq. H82|) for the i < 2 spherical harmonic components of 
this equation, and a regular freezing boundary condition, 
Eq. 179() . for the i > 2 spherical harmonic components. 

Finally, we note that boundary conditions for and 
Zf can also be obtained from the 'T-freezing" condition 
that is often used to determine the shift vector [s^, HI] . 
We found that these F-freezing boundary conditions on 

and zf are not as effective as Eqs. (|5n)l and 
in controlling the growth of instabilities in these "gauge" 
fields. And since the equations for the F-freezing forms 
of these boundary conditions are quite lengthy, we do not 
reproduce them here. 



VII. NUMERICAL METHODS 

In this section we describe briefiy the numerical meth- 
ods used to compute the simulations presented below. 
All our numerical computations are performed using a 



multidomain pseudospectral collocation method. Our 
numerical methods are essentially the same as those 
we applied previously to evolution problems with scalar 
fields HHil, with the MaxweU system and the Ein- 
stein system [H El IH . 



A. Spectral Collocation Method 

The computational domain for the single black-hole 
evolutions described in Sec. IVIlH is a spherical shell ex- 
tending from some rmin (just inside the black-hole event 
horizon) to some maximum value rmax- This domain 
may also be subdivided into one or more spherical-shell 
subdomains. Consider a single one of these subdomains 
that extends from radius r'^[^ to r^j^^. Given a system 
of partial differential equations 

atu"(x,t) =.F"[u(x,i),9^u(x,i)], (83) 

where is a collection of dynamical fields, the solu- 
tion (x, t) on this subdomain is expressed as a time- 
dependent linear combination of N spatial basis functions 



N-l 



.^(x,i)= ^<(O0fc(x). 



(84) 



k=0 



We expand each Cartesian component of each tensor 
in terms of the basis functions Tn{p)Yim{9,^), where 
y«m(^, ¥') are spherical harmonics and Tn{p) are Cheby- 
shev polynomials with 



2r- 



P = 



(85) 



The spherical coordinates (r, 9, and tp) used in these 
spectral expansions are related to the pseudo-Cartesian 
coordinates used to evaluate the components of tensor 
fields {e.g. gxx, Kyz, D^yz, •■■) by the usual transforma- 
tion X = r sin 9 sin ip, y = r sin 9 cos if, and z — r cos 9. In 
these spectral expansions we keep Chebyshev polynomi- 
als Tn{p) with n up to some finite Nr^ and spherical har- 
monics Yim{9, <f) with i up to some finite fmax- For each 
£, we keep all \m\ <£. The values of Nr and ^max deter- 
mine our radial and angular resolution, and we vary these 
in order to perform convergence tests. Spatial derivatives 
are evaluated analytically using the known derivatives of 
the basis functions: 



N-l 



dMx,t) = ^ <(i)9,0fe(x). 



(86) 



k=0 



Associated with the basis functions is a set of Nc col- 
location points Xi. Given spectral coefficients u'^{t), the 
function values at the collocation points u" (x^ , t) are 
computed using Eq. H84() . Conversely, the spectral co- 
efficients are obtained by the inverse transform 



E 

i=0 



WlU^(x^,t)0fe(Xj), 



(87) 
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where Wi are weights specific to the choice of basis func- 
tions and collocation points. Thus it is straightforward 
to transform between the spectral coefficients u^{t) and 
the function values at the collocation points u'^{'K.i,t). 
The partial differential equation, Eq. (|83|1 . is now rewrit- 
ten using Eqs. H84|l - H87|l as a set of ordinary differen- 
tial equations for the function values at the collocation 
points, 

5tu^(x„i) =gn«iv(x,,t)], (88) 

where Qf depends on u'^{xj,t) for all j. We integrate 
this system of ordinary differential equations, Eq. (|88|l . 
in time using a fourth-order Runge-Kutta algorithm. 
Boundary conditions are incorporated into the right side 
of Eq. using the technique of Bj0rhus j^^l- For the 
evolutions reported here the time step is typically chosen 
to be about 1.5 times the distance between the closest 
collocation points in order that the Courant-Friedrichs- 
Lewy stability limit remains satisfied. 

We use no filtering on the radial basis functions, but 
apply a rather complicated filtering rule for the angular 
functions. When evaluating the right side of Eq. (|88|l . 
we set to zero the coefficients of the terms with £ > 
^max — 3 in the tensor spherical-harmonic expansions 
of the dynamical-field components {i.e. the gij, Kij, 
and Dkij components) of this equation. This filtering 
method eliminates a certain type of angular instabil- 
ity that arises because differentiation mixes the various 
spherical-harmonic indices in the spectral expansions of 
the Cartesian components of tensors. 

B. Multiple Subdomains 

Under many circumstances it is advantageous to divide 
the computational domain into multiple subdomains. 
This subdivision allows faster calculations by spreading 
the computational load over multiple processors, ft also 
allows the use of different spectral resolutions in differ- 
ent parts of the computational domain, thus making it 
possible to concentrate computational resources in ar- 
eas where they are most needed. We use the spectral 
collocation method described in Sec. lVff Xl for each com- 
putational subdomain. An additional complication with 
multiple subdomains is that only some of the subdo- 
main boundaries are external, while some are just in- 
ternal boundaries that separate subdomains. Appropri- 
ate boundary conditions must nevertheless be specified 
on each boundary of each subdomain. For the computa- 
tions described here we impose conditions on each subdo- 
main boundary using the Bj0rhus [s^ method, in which 
dtu" [see Eq. (H^ ] is specified for each ingoing charac- 
teristic field {i.e., for each d with V(^a) < 0)- Non-ingoing 
fields do not need and are not given boundary conditions. 
For external boundaries (those without neighboring sub- 
domains) the incoming dtu" are computed according to 
some externally imposed boundary condition, for exam- 
ple our new constraint-preserving or physical boundary 



conditions. For internal boundaries, the values of all in- 
coming dfu" are simply copied from the corresponding 
outgoing dtu" in the neighboring subdomain. At a black- 
hole excision boundary V(^a) > for Q;, so no boundary 
condition is required there on any of the dynamical fields. 

For the simulations described here, we subdivide our 
computational domain into a set of concentric spherical- 
shell subdomains. We choose the constants Nr and £,nax 
(that specify the spectral resolution of each subdomain) 
to have the same values in all the subdomains; this makes 
it easier to achieve load balancing when multiple proces- 
sors are used. We concentrate the numerical resolution 
where it is needed in our solutions by allowing spherical 
shells with varying thickness: for fixed radial resolution 
Nr, thin shells achieve higher resolution than thick shells. 

C. Residual Evaluators 

We use three different computational diagnostics to 
help us evaluate the accuracy and stability of the nu- 
merical simulations described in Sec. IVfffI The first and 
simplest diagnostic just measures the difference between 
a numerical solution and the exact solution of the equa- 
tions. This is possible whenever the exact solution to the 
problem is known, as in our evolutions of unperturbed 
black-hole spacetimes. We measure the deviation of the 
numerical solution from the exact solution quantitatively 
by evaluating the energy norm of the difference between 
the two solutions: 




Here denotes the numerical solution, Uq the exact so- 
lution, and Saf3 the symmetrizer matrix of the hyperbolic 
evolution system. The symmetrizer 5*0.^ is a positive def- 
inite matrix on the space of dynamical fields, which is 
derived for the KST evolution system studied here in 
Refs. [s^, E3- (We evaluate Sai3 using the exact solu- 
tion Uq.) The energy 6E is not dimensionless, and con- 
sequently it is not clear how to interpret whether a given 
value means that the solution is a good approximation 
of the exact solution or not. Therefore we normalize by 
dividing by the total energy of the exact solution: 

^0 = J So,pu'^4^d^x. (90) 

The ratio SEjE^ is therefore a good dimensionless mea- 
sure of the accuracy of our numerical solution. When 
this ratio becomes of order unity, the numerical solution 
bears little resemblance to the exact solution. 

A second, more general, method of measuring the 
accuracy and stability of our numerical solutions is 
to monitor the magnitudes of the constraints, — 
{C,Ci,Ckij,Ckiij}. When the constraints vanish, the solu- 
tions to our first order evolution system, Eqs. (0-©, are 
guaranteed to be solutions of the original Einstein equa- 
tions as well. Thus, to measure how well our numerical 
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solutions solve the original Einstein system we construct 
the following measure of the constraints: 

IICIP = JiC^+C.C+Ck^f"'^ +CfeH,C"*^)V?rf'2;.(91) 

The constraints consist of various derivatives of the dy- 
namical fields . Therefore we can construct a meaning- 
ful dimensionless measure of the constraints by normal- 
izing ||C|| by dividing by a measure of these derivatives: 



\\du\ 



(92) 



We evaluate the quantities that appear in Eq. {i.e., 
the tensor components of u" and the partial derivatives) 
in the global Cartesian-like coordinate system used in our 
code. The ratio ||C||/||9u|| is therefore a dimensionless 
measure of the degree to which our numerical solutions 
satisfy the original Einstein system. When this quantity 
becomes of order unity the constraint violations domi- 
nate and our numerical solutions no longer accurately 
represent solutions to the original Einstein equations. 

At the analytical level, the vanishing of ||C|| is enough 
to ensure that the solutions to our first-order evolution 
system also represent solutions to the original second- 
order Einstein equations. At the numerical level, how- 
ever, there remains a (small we think) possibility that 
our solutions to the discrete first-order evolution system 
manage to keep ||C|| small while somehow failing to sat- 
isfy the original second-order Einstein system. In order 
to detect whether this has occurred, we construct a third 
independent measure of the degree to which our numer- 
ical solutions satisfy the original system. Following the 
ideas of Choptuik '43l| , we use our numerical solution for 
to reconstruct the four-dimensional metric: 

ds^ — g^j^udx^dx^ 

= -N^dt^ -f g,j{dx' -f N'dt){dx^ + N^dt). (93) 

[Letters from the latter part of the Greek alphabet (/z, 
V, etc.) denote four-dimensional spacetime indices in 
our Cartesian-like coordinate system.] We evaluate the 
components of g^i, on the actual time slices used in our 
evolution code, and then evaluate the time derivatives 
dtg^ii^ and d^g^i, using standard centered finite-difference 
expressions for these quantities with seven-point sten- 
cils. Given these time derivatives on a certain time slice, 
we then evaluate the complete set of first and second 
derivatives, d^g^u and da-drg^u, by computing the spa- 
tial derivatives using the spectral methods described in 
Sec. IVII Al Finally, we combine these first and second 
derivatives to determine the four-dimensional Ricci ten- 
sor: 

Rf^iy = -^^g"'' {dcrdrgt_iu^d^,d„gaT-'2^dad{f_,g^)r) 



<JT UJp 
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df_cgawd^grp + '^dag^,wd[rgp\v 

{d^gt,^ - 29(^3^)^) {dpg^r - '^d^grp) ■ (94) 



The Ricci tensor should vanish for any solution of the 
vacuum Einstein equations, so this quantity gives us an 
independent way to measure how well our numerical solu- 
tion solves the original second-order Einstein equations. 

The expression for the Ricci tensor, Eq. 194|) . has thir- 
teen terms (each of which contains contractions): four 
terms proportional to dadrg^^, and nine terms propor- 
tional to dagp,v (when symmetrizations and antisym- 
metrizations are expanded out). We construct a positive 
definite quantity, , to which the size of i?^^ can be 
compared, by taking the square root of the sum of the 
squares of these thirteen terms (for each fi and v). We 
then construct norms of these quantities by integrating 
them over our spacelike slices: 



\R 



\R\ 



RMS I I 2 



Y,{R^.uf^/gd' 



Y.^RT^"f^/gd' 

/J,<Z-' 



(95) 
(96) 



We note that the i?™^ 



and used to compute these 
norms are the quasi-Cartesian components of these quan- 
tities used in our code. We use the ratio | |/| | 
as our third independent measure of the degree to which 
our numerical solutions satisfy the original Einstein equa- 
tions. When both 1 1 i? 1 1 / 1 1 i?"*'" 1 1 and 1 1 C 1 1 / 1 1 9m 1 1 are small , 
we have considerable confidence in the accuracy of our 
solutions even in cases where the exact solution is not 
known. 



VIII. NUMERICAL EVOLUTIONS 

We have studied the efficacy of the boundary condi- 
tions proposed in Secs. lIVtfVll bv performing evolutions of 
various single black-hole spacetimes. These simulations 
evolve either unperturbed Schwarzschild black holes in 
Kerr-Schild coordinates. 



2M 



ds^ :^-dt' + (dt + dry+ dr^ + r'dVL\ (97) 

r 

or fully dynamical black holes that are small perturba- 
tions of the Schwarzschild spacetime (but we solve the 
full nonlinear equations in our code). We express all 
lengths and times associated with these simulations in 
units of the bare black-hole mass Af , even for the per- 
turbed solutions in which the ADM mass exceeds M . 
The computational domains for these evolutions con- 
sist of one or more concentric spherical shells that cover 
the space from rmin = 1.9Af (just inside the black-hole 
event horizon) to some maximum value rmax- We have 
run evolutions using several different Tmax in the range 
6.9Af < r^ax < 41.9Af. 

The initial data for these evolutions are prepared by 
applying an odd-parity outgoing quadrupole-wave per- 
turbation with amplitude 4 x 10~^ to the Kerr-Schild 
spatial metric and its time derivative. This outgoing 
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FIG. 1: Initial data for the perturbed black-hole evolutions. 
Shown is the Newman-Penrose component of the Weyl tensor 
1^*4 1, which measures the outgoing gravitational wave flux. 



wave pulse is constructed by Teukolsky's method |4J| 
with the generating function G{r) that determines the 
radial profile of the wave: G{r) — Ae^^'^~^°'> , where 
A = 4x10^'^, To = 5M and w = 1.5M. Then we solve nu- 
merically the full non-linear conformal "thin-sandwich" 
form of the initial value equations to obtain constraint- 
satisfying perturbed black-hole initial data (see Pfeif- 
fer, et al. ^3). In the resulting initial data set, the 
ADM energy exceeds the apparent horizon mass by about 
10~^M. Figure ^ illustrates the initial data for these 
perturbed black-hole evolutions. The quantity shown in 
Fig. n is a measure of the outgoing gravitational wave 

flux, ^Uf+Ul+g^^gi^ /2, defined in Eq. GH)- This quan- 
tity is equivalent to the Newman-Penrose component of 
the Weyl tensor |^'4|. The outer radius of the compu- 
tational domain for the initial data shown in Fig. ^ is 
r-max = 11.9M. 

For simplicity we use fixed gauge evolutions: the lapse 
density Q and the shift iV* are set to their initial Kerr- 
Schild values, Eq. (|97|l . for all times. For the perturbed 
initial data, we also set Q and to the unperturbed 
analytical Kerr-Schild values, ignoring the solutions for 
these fields given by the conformal "thin-sandwich" ini- 
tial value equations. Although more sophisticated gauge 
conditions should result in better long-term evolutions, 
these simple fixed-gauge conditions are adequate for the 
tests of the new constraint preserving boundary condi- 
tions, which are of primary interest here. 

The KST evolution system has a number of freely 
specifiable parameters. We set the KST parameters to 
the values 70 = 0.5, 71 = —12, 72 = —1, 73 = 0.16, 
and 74 = —0.96 for all the evolutions discussed here. 
This choice of parameters is one of the special cases 
that leave the parameter q of Eq. (|A14ll ill-defined. In 
this case the choice of q is arbitrary, and we set q — 1 
in all of the evolutions discussed here (since this yields 
better performance than q = 0). For these parameter 



choices, both the fundamental evolution equations and 
the constraint evolution system are symmetric hyper- 
bolic, and all characteristic speeds (relative to hypersur- 
face orthogonal observers) are either or ±1; in particu- 



lar, VI 



1. The values of the KST parameters 



used here are the same as those used in previous stud- 
ies 0,^3. However, the current evolution equations are 
rather different from those used in Refs. 40, 41J, in which 
the evolution equations were modified by a kinematical 
change of variables and by the addition of terms propor- 
tional to the constraint Ckij- Here we perform no such 
modification; that is, the kinematical parameters defined 
in Refs. [2^ liol l4l| are chosen here to have the values 
k = l,z = a = h = c = d = e^Q. The values of the KST 
parameters in Refs. [40ll4l| were chosen to minimize the 
growth rate of the constraints for evolutions of a single 
black hole in Painleve-GuUstrand coordinates. Because 
we use a rather different set of evolution equations, and 
because we evolve black holes in Kerr-Schild coordinates, 
we do not expect that these parameter values are optimal 
for our current evolutions. 



A. Unperturbed Black Holes 

In this subsection we describe three numerical tests 
involving 3D evolutions of unperturbed Schwarzschild 
black holes. These evolutions use the Schwarzschild spa- 
tial metric and extrinsic curvature in Kerr-Schild coordi- 
nates, Eq. JHTJ), as initial data. The first and second tests 
explore the evolutions of these black holes using simple 
"freezing" boundary conditions. Freezing boundary con- 
ditions set the time derivatives of each of the incoming 
characteristic fields to zero at the boundaries, i.e., 



dtu^ = 0, 



(98) 



for all a with V(^a) < [where dtu" is defined in Eq. (|49l) ]. 
These boundary conditions are known to make the evo- 
lution equations mathematically well-posed, but they 
are also known to do a poor job of preventing the in- 
fiux of constraint violations into the computational do- 
main for fully dynamical time-dependent solutions 1 , 12|| . 
Freezing boundary conditions, however, have been used 
rather successfully for time independent solutions such 
as the unperturbed black-hole spacetimes considered in 
this subsection poll4ll| . 

Our first test evolves an unperturbed Schwarzschild 
black hole on a computational domain extending from 



1.9A/ to Tn 



41. 9M. This domain is sub- 



divided into eight subdomains, each of width 5Af , each 
using the same angular resolution ^max = 21, and each 
having the same radial resolution Nr = 17, 21 or 26. 
This angular resolution is much higher than necessary 
to resolve this spherically symmetric spacetime, but we 
wanted to verify the stability of our code even for large 
values of ^max- Figure 121 shows two measures of the errors 
in these evolutions: ||C||/||9m|| and (5i?/i?o. These results 
show that our error measures converge toward zero as we 
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FIG. 2: Constraint violations and energy norms for unper- 
turbed black-hole evolutions using freezing boundary condi- 
tions. The outer boundary is at rmax = 41. 9M; and three 
different radial resolutions, A^'r — 17, 21, and 26, show nu- 
merical convergence. 



increase the radial resolution Nr. These results also show 
that our computational methods are capable of evolv- 
ing unperturbed black-hole solutions for hundreds (if not 
thousands) of M, which is consistent with previous re- 
sults using a similar evolution system |43, 113 • 

As a second test, we explore the effects of changing the 
location of the outer boundary of the computational do- 
main, Tmax, for thcsc Unperturbed black- hole evolutions 
with freezing boundary conditions. Figure 13 shows the 
constraint violations ||C||/||9u|| and the energy norms 
5E / Eq for evolutions using several different outer radii: 
r^^^ = 6.9M, U.9M, 16.9M, and 41.9M. All of these 
evolutions use the same ^max = 21, and Nr ~ 26 in each 
subdomain. Each computational subdomain has width 
5M, so the number of subdomains is adjusted to achieve 
the desired rmax- Figure |3| shows that these evolutions 
have an instability that causes ||C||/||9w|| and 6E/Eq to 
grow exponentially. This instability becomes weaker as 
''max increases; we find (by measuring the growth rate for 
r^^^ = 6.9M, 11.9M, 16.9Af, 21.9M, and 31.9M) that 
the growth rate is M/r « g-''max/i3M^ ^^lig appears to 
be a constraint-violating instability that is influenced by 
the location of the outer boundary of the computational 
domain. Figure 01 demonstrates for the r,„ax = 11. 9A/ 
case that these unstable evolutions are numerically con- 
vergent, by showing that ||C||/||(9u|| and SE/Eq approach 
an exponentially growing solution as Nr increases. Fig- 
ure 21 also shows that the growth rate of the instability is 
the same for all values of Nr in this r,„ax = 11.9Af case. 
(Numerical convergence and independence of the growth 
rate with Nr is also observed for the other values of rmax 



FIG. 3: Constraint violations and energy norms for unper- 
turbed black-hole evolutions with a range outer boundary 
radii, nnax = 6.9M, 11. 9M, 16. 9M, 41.9A'/. These runs aU 
use Nr — 26 in each subdomain, but use different numbers of 
subdomains to achieve different rmax. They all use freezing 
boundary conditions, and £max = 21. 



that we tested.) We have also verified that the growth 
rate of this instability does not depend on the angular 
resolution £max for 7 < ^max < 21. These checks suggest 
that this constraint violating instability is a solution of 
the continuum differential equations, and is not primar- 
ily an artifact of the discrete numerical representation of 
the equations used here. 

For the third numerical test we use our new bound- 
ary conditions to evolve the same unperturbed black-hole 
initial data on the same computational domains used in 
the first two tests. These new boundary conditions in- 
clude the new constraint preserving boundary conditions, 
Eqs. (O-linnil with ^5 = 0.75 and /xe = -0.5, the new 
physical boundary condition, Eq. lf75)l . and the new gauge 
boundary conditions, Eqs. (|80|l and H82|) . Figure|31depicts 
||C||/||9m|| and SE/Eq for evolutions of the unperturbed 
black-hole initial data on a computational domain with 
rmax — 41.9Af. (Except for boundary conditions, the 
evolutions in Fig. [3 are identical to those of Fig. |21) The 
evolutions of Fig. |31 show that both the constraints and 
the energy norms of these solutions decrease toward zero 
as Nr increases. 

At late times, the highest resolution curve plotted in 
Fig. |S1 shows the beginning of an exponentially-growing 
mode. We call this a gauge mode, because SE/Eq be- 
gins to grow in these solutions without a corresponding 
growth in the constraints ||C||/||9m||. This gauge mode 
grows more rapidly if we move the outer boundary in- 
ward, as can be seen in Fig.|^ which shows the analogous 
evolutions with outer boundary at rmax = 11.9Af. The 
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FIG. 4: Constraint violations and energy norms for an un- 
perturbed black-hole evolution using freezing boundary con- 
ditions. The outer boundary is at rmax = 11. 9M; and three 
different radial resolutions, Nr = 17, 21, and 26, show numer- 
ical convergence to an exponentially growing solution. 



FIG. 5: Constraint violations and energy norms for unper- 
turbed black-hole evolutions using our new boundary condi- 
tions. The outer boundary is at rmax = 41. 9M; and three 
different radial resolutions, Nr — 17, 21, and 26, show nu- 
merical convergence. 



constraints remain roughly constant in the evolutions of 
Fig. El for a considerable amount of time after SE/Eq 
begins to grow rapidly. Comparing the results of Fig. 01 
with those of Fig.jSl we see that our constraint-preserving 
boundary conditions do improve the constraint violations 
in these evolutions, but only slightly. However, this im- 
provement comes at the expense of introducing a new 
gauge-mode instability. At very late times we see that 
the highest resolution evolution also shows signs of a 
constraint violating instability, although it grows more 
slowly than the gauge mode. The growth rate of the con- 
straints seen in the highest resolution evolution of Fig.O 
is the same as the growth rate of the instability seen in 
Fig. ^ to within about 5%. This suggests that both of 
these constraint violating instabilities might be caused 
primarily by bulk rather than boundary constraint vio- 
lations. 

The gauge-mode instability seen in Figs. [3 and |S1 is 
dominated by its ^ = 3 spherical harmonic component. 
Recall that we impose the gauge-fixing boundary condi- 
tion, Eq. H82|l . by filtering out everything but the £ < 2 
components of dtU^~ . If we change this filtering to im- 
pose Eq. (|82|l on only the £ < 1 components of dfU^^ , 
the gauge-mode instability is dominated by its £ = 2 
spherical harmonic component, and grows more rapidly 
than in Figs. [SI and |S1 Thus the gauge boundary con- 
dition, Eq. (|S^ . does improve the stability of the gauge 
mode at least to some degree. However, if we impose 
Eq. lIS2Il without any filtering, the instability grows much 
faster than the rate seen in Figs.ElandjSl and this growth 
rate increases as fmax increases, so the evolution becomes 



non-convergent. Recall that we use the simplest possible 
gauge conditions: we set the lapse density Q and the shift 
A^* to their analytic values throughout the evolution. A 
more sophisticated treatment of gauge conditions (not 
just at the boundaries, but throughout the volume) is 
probably needed to control these unstable gauge modes. 



B. Perturbed Black Holes 

In this subsection we describe 3D numerical evolu- 
tions of perturbed black-hole spacetimes. The initial 
data for these evolutions are discussed at the beginning 
of Sec. IVIIII We evolve these initial data using the 
same non-linear equations and evolution methods used 
in Sec. IVIII Al for the unperturbed cases. These per- 
turbed black-hole initial data include a short-wavelength 
gravitational- wave packet, so more radial collocation 
points are required to achieve an accuracy comparable 
to that of the unperturbed black-hole cases. 

Figure [71 shows the constraint error ||C||/||9m|| for evo- 
lutions of these perturbed black-hole initial data. These 
simulations are performed on a computational domain 
with four concentric subdomains, each of width 5M and 
^max = 11, and with outer boundary at rmax = 21.9. The 
dashed curve in Fig.[71shows the results of using the sim- 
ple freezing boundary conditions, Eq. a large con- 
straint violation is generated at t « 20M when the wave 
pulse passes through the outer boundary of the compu- 
tational domain. The size of this constraint violation 
is comparable to the amplitude of the wave, and does 
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FIG. 6: Constraint violations and energy norms for an un- 
perturbed black-hole evolution using our new boundary con- 
ditions. The outer boundary is at rmax = 11.9M; and three 
different radial resolutions, Nr = 17, 21, and 26, show numer- 
ical convergence to an exponentially growing solution. 



FIG. 7: Constraint violations for evolutions of perturbed 
black holes using freezing (dashed curve) and our new bound- 
ary conditions (solid curves). Radial resolutions A'^ — 21, 
31, 41, and 51 are shown for the new boundary conditions, 
while only A''^ — 51 is shown for the freezing case. The outer 
boundary is at rmax ~ 21.9. 



not converge away with increased resolution. Only the 
highest radial resolution is shown, because the lower res- 
olution curves are almost identical. Therefore, these sim- 
ulations with freezing boundary conditions fail and can 
not be used to model the physical gravitational waves in 
these solutions after the time t ~ 20M with any accuracy. 

Also plotted in Fig. [7| are evolutions of the same per- 
turbed black-hole initial data on the same computa- 
tional domain, but now using our new boundary con- 
ditions: the new constraint preserving boundary condi- 
tions, Eqs. with ^5 = 0.75 and = -0.5, 
the new physical boundary condition, Eq. lf75|l . and the 
new gauge boundary conditions, Eqs. H80|l and H82|l . In 
these cases the constraint violation generated as the wave 
passes through the outer boundary converges away with 
increasing radial resolution. The constraint violation 
is eight orders of magnitude smaller than with freezing 
boundary conditions in the highest resolution case. Fig- 
ure|S|shows an independent measure of how well these nu- 
merical evolutions satisfy the original Einstein equations 
by plotting the average value of the four-dimensional 
Ricci tensor jggg ^j. ^j^g ^^^^ 

simulations shown in Fig. [7| Figure |S1 confirms the re- 
sults seen in Fig. [7| the Einstein equations are violated 
when the wave hits the boundary with freezing bound- 
ary conditions, while the new boundary conditions are 
quite effective in reducing this violation by many orders 
of magnitude. 

Two interesting features are seen in the evolutions in 
Figs. [71 and |S1 which use the new boundary conditions: 



Freezing BC 




t/M 

FIG. 8: Four-dimensional Ricci-tensor residual of Eq. H95^ for 
the perturbed black-hole evolutions shown in Fig. |7| 



first, there is an unstable exponential growth in the high- 
est resolution runs starting at about t « 200M, and 
second, there is some complicated structure starting at 
about t « 20M which is reduced and finally vanishes at 
the highest radial resolution. Consider first the unsta- 
ble exponential growth. Figure |51 shows this exponential 
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FIG. 9: Constraint violations in evolutions of a perturbed 
black hole. The dashed curve corresponds to freezing bound- 
ary conditions, and the solid curves correspond to the new 
boundary conditions for different rmax- 



growth more clearly by displaying the constraint viola- 
tions for a range Tmax- The simulations in Fig. |51 have 
iVr = 51 and ^max = 11 in each subdomain, while rmax is 
varied by adjusting the number of subdomains. However, 
in order to resolve the small incoming waves reflected 
off the outer boundary, we need higher resolution near 
the outer boundary. The subdomains are made narrower 
therefore for larger r, as described below. Figurel^shows 
that the exponential growth rate of this constraint violat- 
ing instability decreases as rmax increases. The growth 
rate of this constraint- violating instability is about 1.7 
times larger than the growth rate seen in the unperturbed 
black-hole evolutions of Fig. |2| for each rmax- The simi- 
larity in size and dependence on the location of the outer 
boundary suggests that the instabilities in Figs. |3| and El 
may have the same basic cause. If so then this cause is 
probably bulk generated constraint violations, since the 
evolutions in Fig. |51 use constraint preserving boundary 
conditions. 

The second interesting feature seen in Figs. [71and|Hlis 
the complicated structure of the curves starting at about 
t w 20M, when the gravitational wave passes though the 
outer boundary. This structure seems to be caused by 
the (very weak) reflection of the gravitational wave pulse 
back into very short wavelength ingoing waves. Consider 
an outgoing wave with coordinate velocity near one (the 
speed of light) as it approaches the outer boundary. Since 
the non-linear evolution equations and our boundary con- 
ditions couple the various characteristic fields, this wave 
will be mixed with and partially reflected as an ingoing 
wave which propagates at the much smaller shift speed 
N'^rii <C 1, even in the continuum limit. If the original 



outgoing wave has wavelength A, then the reflected in- 
coming wave that propagates at the shift speed will have 
the much smaller wavelength XN^Ui <C A. As the outer 
boundary radius is increased, the amplitude of these re- 
flected waves is decreased. But their wavelength also 
decreases (because the shift decreases as r increases) 
and therefore these reflected waves become more difficult 
to resolve numerically for larger rmax- 

One approach is to ignore these very small amplitude 
non-physical reflected waves and not even attempt to re- 
solve them. However, these reflected waves contribute 
(slightly) to the constraint quantities, so leaving them 
unresolved would introduce constraint violations that are 
roughly the size of the reflected waves. If this constraint 
violation is smaller than numerical truncation error in the 
remainder of the domain, then no harm is done by ignor- 
ing the reflected waves. However, for the pseudospec- 
tral simulations presented here, the truncation error is 
so small elsewhere that the contributions of the reflected 
waves can dominate, obscuring our convergence tests. 
Therefore, we choose to resolve the reflected waves in 
the tests presented here. 

Another approach would be to attempt to eliminate 
the problem completely by making a smarter choice for 
the shift N'^ . The radial component of the shift could be 
made to approach a constant value rather than falling 
to zero as r increases. This would limit the amount 
of blue-shift the reflected waves could experience. Al- 
ternatively the radial component of the shift could be 
made to pass through zero and switch sign ^18J, thus 
eliminating the shift-speed incoming waves completely. 
These changes would either limit or entirely eliminate 
the problem, but possibly at the expense of introduc- 
ing other gauge-related difflculties. Since the choice of 
gauge is not the main subject of this paper, we decided to 
deal with the reflected-wave problem here by increasing 
the resolution in the subdomains near the outer bound- 
ary. For instance, the rmax = 26.9M curve in Fig. ^ 
was produced using eight subdomains with boundaries 
at r = 1.9M, 6.9M, 11.9M, 16.9M, 19.4M, 21.9M, 
23.567M, 25.233M, and 26.9M, each using Nr = 51 col- 
location points. 



C. Angular Instability 

In addition to the constraint violating instability and 
the gauge mode instability discussed in Sees. I Vlll Al and 
IVIII Bl there are also signs of a very weak instability 
that primarily affects the highest angular modes of the 
evolved flelds. This angular instability has a growth rate 
that increases with increasing ^max, and so this insta- 
bility appears to be non-convergent. This instability is 
not evident in any of the flgures shown so far, and is 
negligible on the time scales of interest here except for 
simulations using very large ^max and very small rmax • In 
order to see this instability clearly, we look at a quantity 
that is linear in the dynamical fields and vanishes unless 



16 




FIG. 10; Highest unfiltered components of SKf^^, defined in 
Eq. lHUTt . for evolut ions of unperturbed black holes. These 
evolutions use a single subdomain with Nr = 26 and rmax ~ 
6.9M. 



FIG. 11: SRf'^^ for unperturbed black-hole evolutions on 
domains with different rmax. This quantity is evaluated on the 
outermost domain which has width 5M, A''^ = 26, and ^max ~ 
21. These evolutions use our new boundary conditions. 



the instability is present. To this end we define 

6K = g'^K,j-K,j), (99) 

where g^^ and Kij are exact solutions for the three-metric 
and extrinsic curvature, and Kij is the numerical extrin- 
sic curvature from our simulations. To see the angular 
instability we project SK onto the spherical harmonic 
basis at each radial collocation point rp, 

SKpirn^ J Yi^{e,^)SK{rp, e,^) sin edOd^, (100) 
and then average these over Vp and m by forming 

{5Kr''f^N-\2£+ir' ^ {6Kp,,^r. (101) 

p. \m\<i 

Figure [TUl shows the highest unfiltered components of 
SK^^^ for evolutions of unperturbed black holes. These 
evolutions use our new boundary conditions, and are 
performed on a single computational subdomain with 
Nr — 26 and rmax = 6.9M. Each curve in Fig. ^|is gen- 
erated from a run with a different imax, and the SKf^^ 
plotted is the one with £ = fmax — 4, the largest i that 
is untouched by our angular filtering procedure. Each £ 
component grows exponentially at a rate that increases 
with £. (The extremely rapid growth in the £ = 10, 12, 
and 15 components that occurs just before the simula- 
tion crashes is presumably due to nonlinear coupling that 
becomes important only when the gauge mode becomes 
large.) For all cases plotted in Fig.^|except £max = 21, 
the simulation crashes because of the (convergent) gauge 



mode instability described earlier. For these cases the an- 
gular instability is orders of magnitude smaller than the 
unstable gauge mode for the entire duration of the sim- 
ulation. Only for fmax = 21 does the angular instability 
dominate before the end of the simulation. 

We have explored the behavior of this angular insta- 
bility in two ways. First we verified that the growth rate 
for a given f-component of SKf^^ is independent of the 
^max used to computc it. We did this by comparing the 
curves in Fig. 1 101 with graphs of the same quantities com- 
puted from a single run with ^max = 21. The resulting 
plot looks almost identical to Fig.^1 except at very late 
times when the simulations begin to crash. And second, 
we explored the growth rates of the angular instability for 
fixed ^ as a function of rmax- Figure HI] shows SKf^^ for 
runs with different rmax- We increased r^ax in these runs 
by adding more subdomains of width 5M, each having 
Nr — 26 and £max = 21. Because the angular instability 
is largest near the outer boundary, we compute SK^^^ 
only in the outermost subdomains for these plots. 

The curves in Fig.^Jshow that the angular instability 
becomes weaker and weaker as we increase the size of the 
computational domain. To study this behavior quantita- 
tively, we plot the exponential growth rates of SKff^^ as 
a function of rmax; these are shown as circles in Fig. 1121 
The error bars in the inset correspond to the variation 
in slopes that we extract from data segments of different 
length, but these error bars are not shown in the main 
plot in Fig. ^] because they are the same size as the plot 
symbols. Rather than falling off gradually like l/rmax 
as might have been expected, the growth rate in Fig. 1121 
appears to go to zero at a finite value of rmax- This finite 
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FIG. 12: Exponential growth rate of SKf^^ (circles) and 
5A'i2^^ (triangles) for different rmax in evolutions of unper- 
turbed black holes with our new boundary conditions. The 
inset shows an enlargement with the best linear fits through 
the points 16.9M < r^ax < 26.9M. 



value of Tinax depends on i; to show this dependence, 
we also plot in Figure ^1 the exponential growth rates 
of 6X^2^^ (shown as triangles in the plot) for the same 
runs. The best linear fits through the ^ = 17 and £ — 12 
growth rates have the same slope, but the growth rate 
of the £ = 12 mode goes to zero at a smaller r^ax- Fig- 
ure El suggests that for a given £max there are no angular 
instabilities at all when r,„ax is sufficiently large. No an- 
gular instability has been detected in any of our runs with 
^max <21,t< 300M, and r^ax > 27M. Our results also 
suggest that for a given rmax, an angular instability with 
arbitrarily large growth rate could be found by making 
^max sufficiently large. 

We see no angular instability at all when we use freez- 
ing boundary conditions. Furthermore, the angular in- 
stability shown in Figs. I10H12I is present whether or not 
we use the physical boundary condition, Eq. l(7S|) . or 
the gauge boundary conditions, Eqs. llSn|) and [Al- 
though as noted before, an angular instability dominates 
when Eq. H82() is imposed on all the spherical harmonic 
components of U^^ .] It is unclear whether this angular 
instability is due to our numerical method or whether 
the new constraint-preserving boundary conditions yield 
an ill-posed initial-boundary-value problem at the con- 
tinuum level. We note however that for the resolutions 
and time-scales of interest here, this instability remains 
small and can be controlled quite effectively by modestly 
increasing the radius of the outer boundary. And the nu- 
merical evolutions, such as those in Figs. [THOl produced 
by these methods do appear to be accurate solutions of 
the Einstein equations: both the constraints ||C||/||9u|| 



and the four-dimensional Ricci tensor | |i?| |/| | can 
be made arbitrarily small. Therefore from a practical 
point of view it may not matter whether these bound- 
ary conditions are formally well-posed, or that our com- 
putational methods contain a very mild non-convergent 
angular instability. 



IX. DISCUSSION 



This paper constructs new boundary conditions for the 
KST form of the Einstein evolution system that are de- 
signed to prevent the influx of constraint violations and 
physical gravitational waves into the computational do- 
main. From a mathematical point of view, these bound- 
ary conditions are Neumann-like (in the sense that they 
place conditions on the normal derivatives of the in- 
coming characteristic fields). Boundary conditions of 
this type have not been studied as comprehensively as 
the simpler "maximally-dissipative" boundary conditions 
(which are Dirichlet-like in that they place conditions on 
the values of incoming characteristic fields themselves). 
Rigorous mathematical well-posedness theorems do not 
yet exist for hyperbolic evolution systems with these 
Neumann-like boundary conditions. So additional math- 
ematical analysis of these conditions is urgently needed 
to determine whether they are well-posed both at the 
continuum differential equation level and the discrete nu- 
merical level. Should this analysis reveal that these con- 
ditions are ill-posed, then alternate ways of preventing 
the influx of constraint violations and physical gravita- 
tional waves in these systems would be needed even more 
urgently. 

Our numerical tests of these new constraint preserv- 
ing and physical boundary conditions show them to be 
quite effective: constraint violations can be reduced to 
roundoff- level errors in dynamical black-hole evolutions. 
Nevertheless, several weak instabilities did appear in our 
numerical results, and additional work is needed to sort 
out exactly what their causes are and what methods 
can be used to control them. Are the constraint vio- 
lating instabilities seen in the evolutions of Figs. 13 and 
1^1 really caused by bulk constraint violating terms in the 
constraint evolution equations (as we surmise), and can 
optimal constraint projection methods such as those de- 
veloped for the scalar field system be used to control 
them? Can the gauge instabilities seen in Fig.|Hlbe con- 
trolled by the introduction of dynamical evolution equa- 
tions for the lapse and shift? Are the non-convergent 
angular instabilities seen in Figs. 1101 and 1111 a symptom 
of ill-posedness of these boundary conditions at the con- 
tinuum or the discrete numerical level? And what can 
be done to cure these problems? 
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APPENDIX A: KST CHARACTERISTIC FIELDS 

Explicit expressions are given here for the characteris- 
tic fields of the KST form of the Einstein evolution sys- 
tem. These characteristic fields consist of the collection 
= {g,„Z\ Zf, Zf, Zf, Z^, Zl^, C/f±, [/3±, 

Uf^}, and can be written in terms of the principal evo- 
lution fields, Eq. Q: 

Z^ = 73n*A^ -2(l + 74)n*Df, 



(Al) 
(A2) 

Zt = iPhD]-2PhD]-4.Phn^n^Djku (A3) 
zf = +48v^n'pJ^n'=Ajfc + 274(5 - 972)P^»i?j 

+3(1 - 372 - 470) (4 - 13){P'^D] + Phn'^n'Djki) 
-2(6 + 74)(5 - 9^2)Phn''n'Djki, (A4) 

5 _ f Tja Tjb ua-b\^kT^ (A5) 

(A6) 

iji^ = ±[l + 2vl + (l + 2ji)q]n'Dl 



76 _ pcafc r) 



1 

2 ■' 



P.,P-')n''D,tk, 



2± 



T [1 - 372 + (1 + 271 + 72)9] Tl^i^f, 
= ±2v2n''P'^K,k + {l + 2-fo)P'^D], 



-(1 - 12)P'^D'' + (270 - J2)P'^n''n'D,kl 



u 



3± 



= ±(1 + 27i)n*D,^ T (1 + 271 + l2)n'D^ 
+V3P'^K,,, 



u: 



4± 



'-P^,P 



ab\ 



Kab ± n'^Dkab T (1 + l'2)n''D^ab}k 



(AT) 
(AS) 
(A9) 

(AlO) 



In these expressions the two distinct traces of Dkij are 
written as 



Dl = P^^D., 

p)2 ^ pjkp) 



ijk, 
kij, 



where 



P'^ = -n'n\ 



(All) 
(A12) 

(A13) 



is the projection orthogonal to n;. The quantities fi, V2 
and 173 are defined in Eqs. (| ^ - 110|l . and q is given by 



9 = ..2 \.2 ■ (A14) 



Finally, the projection tensor P^^^ is defined by 



Q 1 

pcab pc p(a pfc) p pc jjab , ^ pc(a pb) p 

^ki] — ^ j ~ ^^y^ + 2 



_lpabpc^_p^^^ _ pc{apb)^^p 



(i^])k- 



(A15) 



These expressions are the completely general forms of the 
characteristic fields for the KST system. They reduce to 
the expressions for the restricted case vf — V2 — — 1 
published previously in Ref. |29| . 

The characteristic fields U^^, Uf^ and U^^ have char- 
acteristic speeds ±wi, ±z;2, and zLvs respectively (relative 
to the hypersurface normal observers); the fields Ufj^ 
have speeds ±1; the fields {Z^, Zf, Zf, Zf, Zf^, Z^^^} all 
have characteristic speed zero. The characteristic fields 
are linearly independent (so the KST evolution system 
is strongly hyperbolic) if vi ^ 0, V2 ^ 0, V3 ^ 0, and 
7^ '^z- In the case when v\ — v^, the characteristic 
fields given in Eq. (jA7|) are not defined because the 
quantity q given in Eq. (|A14I) is not defined. We find 
that there are nevertheless a complete set of characteris- 
tic fields in the tji = 113 7^ case so long as 



(A16) 



In this case, the characteristic fields C/^^ are just given 
by the expression in Eq. (|A7p with g = 0. Any other 
constant value of q is also acceptable, with the result be- 
ing a redefinition of by the addition of q times l]^^ 
(which has the same characteristic speed in this case). 
The choice g = is probably the simplest, but other 
choices might be desirable under some circumstances. 
For example, for symmetric hyperbolic systems it might 
be better to make the eigenvectors mutually orthogonal 
in terms of the symmetrizer metric. The symmetric hy- 
perbolicity of the KST system is discussed in Appendix 
B of Ref. 113. 

It is also useful to have explicit expressions for the 
inverse transformation, = e"pU^, that expresses the 
principal evolution fields in terms of the characteristic 
fields. These inverse transformations for the general KST 
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form of the Einstein evolution system are given by 



APPENDIX B: HYPERBOLICITY OF 
CONSTRAINT EVOLUTION 



4.V3 

■c ^a^b ] 



(A17) 



Dkij = nuninjn''n°-n Dcab + niUjP'^kn'^n Dcab 
+ (1 + 72)nfe4,) + 2Zf(,n,) + Zl^ 



1 



(A18) 



In this Appendix we evaluate the hyperbolicity of the 
KST constraint evolution system, Eq. (fT^ - ((TH|l . This 
system can be written in the more abstract form 



dtc^ 



t bC 



(Bl) 



In Sec. mil we demonstrated the strong hyperbolicity of 
the KST constraint evolution evolution system by show- 
ing that the characteristic matrices rikA''^ b have a com- 
plete set of eigenvectors. Here we determine the con- 
ditions under which this constraint evolution system is 
also symmetric hyperbolic, i.e., that there exists a sym- 
metric positive-definite Sab (the symmetrizer) with the 
property that it symmetrizes the characteristic matrices: 
A\^^SacA^^b = A%^. 

The most general symmetrizer for the constraint evo- 
lution system can be written conveniently by defining the 
following quantities associated with the constraint CkUj ■ 



The terms involving Dkij on the right sides of Eq. (|A18|) 
are given by 



n^n^u'D^ab = + 271 + 72) + 72(2 + 371 



2vtvi 



77I+ _ Tjl- 

([/3+ _ u-'-) + , , , (A19) 



8vi 



(A20) 



P ^nn D,,,-~^-^Z,+^^(U, + ) 
-[3(l-72)(4-73)-74(5-972)]^ + ^, 



(A21) 



.=i..^^i±p(^--^-)-(l + 27.+72)-., 



(A22) 

(t/^+-C/3-)- (1 + 271):^, (A23) 

^^3 

P-,Dl = (1 - 72) (73 + 274) II - (2 - 372 + 270) II 



2 ^ J3_ 



8vi 

P'.Dl = [73(1 + 270) + 74(2 - 72 + 670)] || 



'4vi 
(A24) 



-(373 -f 274) 



..2 (4-372 + 1470)^ 



2 

(A25) 



'Dkij — £k°'^Cabij, (B2) 

where ekij is the spatial volume element. We also define 



jiki 



(B3) 
(B4) 



+Vlg,,-3Vlg^^,]. (B5) 

Then the most general symmetrizer constructed from the 
metric gij is: 

dS^ = SABdC^dc^ = 

XidC^ + X29"'dCdCa + X35"5^'5'"rf^(fe.,)d^(cab) 

+X5g'''dVldVl + xeg'^'dVfdVl + 2x7g'''dV]dVl. 

(B6) 

The conditions needed to ensure that this symmetrizer 
is positive definite are Xq > fo^' = 1,...,6 and 
X5X6 > X7- There are no cross terms between the mo- 
mentum constraint Ci and the V} or constraints, be- 
cause the resulting terms would have the wrong parity. 
Using the principal parts of the KST constraint evolution 
system given in Eqs. (|15|l ~ p8() . we find that the condi- 
tions needed for symmetric hyperbolicity are 

= Xi(l-^73 +74)-X2(l + 27i), (B7) 
= 2x5(73 + 374) + X27o + X7(274 - 73), (B8) 
= ^X2 -X6(73 - 274)+X7(273 + 674), (B9) 



= ^X272 +X473- 



(BIO) 



20 



The problem now is to determine when these sym- 
metrization conditions can be satisfied. For simphcity 
we focus attention here on the case where all of the char- 
acteristic speeds of the principal evolution system have 
the physical values: 0, ±1. These conditions on the char- 
acteristic speeds implies the following conditions on the 
parameters: 



70 

73 
74 



1 

2' 



-8 



472 -h(5 + 372)(H-27i)' 
1-72- (l + 27i)(5 + 372) 
472 + (5 + 372)(l + 27i) 



(Bll) 
(B12) 

(B13) 



The analysis in Lindblom and Scheel |4C| shows that the 
principal evolution system is symmetric hyperbolic in this 
case for all values of the parameters 71 and 72 that satisfy 
the following inequalities: 



-3 < 72 < 0, 

472 + (5 + 372)(l + 27i) ^0, 



(B14) 
(B15) 



Substituting the conditions Eqs. ljBll|l - ljB13|l into the 
symmetry conditions, Eqs. ljB7|l - ljB10|l . we find the fol- 
lowing conditions on the Xa- 



X2 



Xi 



X5 



X6 



Xi(5 + 372) 



(1 + 271) [472 + (5- 
Xi72(5 + 372) 
16(1-1-271) ' 

Xi(5 + 372) -8x7(1 



■372)(l + 27i)]' 



(B16) 



(B17) 

-f27i)[272 -)- 71(5 -1-372)] 



8(5 + 372)(2 + 77i+672) 
-372)[xi- 8x7(2 + 771 + 67?)] 



(5 



8(1 + 27i)[272 -f 71(5 + 372)] 



(B18) 
(B19) 



The constraint evolution system is symmetric hyperbolic 
iff Xa > for a = {1, 2, 3, 4, 5, 6} and xsXe > Xr- K the 
principal evolution system is symmetric hyperbolic then 
5-1-372 > [from Eq. (|B14|) ]. thus we see from Eq. (|B17|I 
that we must have 



1 -I- 271 < 0, 



(B20) 



or equivalently 71 < — ^ if the constraint evolution sys- 
tem is to be symmetric hyperbolic as well. This condition 
guarantees that 472 + (5 + 372) (1 + 271) < and so the 
second inequality in Eq. (|B15I) is automatically satisfied 
in this case. Consequently Eq. ljB16|l places no further 
restrictions on the parameters 71 and 72. 

To analyze the inequalities on X5 and xe implied by 
Eqs. (|B18|) and l|B19|l . we restrict our attention to the 
simple case where X7 can be set to zero. In this case 
Eqs. (|B18() and l|B19|l simplify considerably: 



X5 
X6 



Xi 



8(27i-hl)(37i+2)' 

(5 + 372)xi 

8(1 + 27i)[272 + 71(5 + 372)] ■ 



(B21) 

(B22) 



Equation l|B22|l guarantees that the coefficient xe is pos- 
itive without any additional restrictions, while Eq. (|B21|) 

shows that the parameter 71 must be further restricted 
by the inequality 



71 < -3, 



(B23) 



in order to ensure that X5 is positive. This argument 
shows then that we can choose the parameters Xa > 
for a = {1, 2, 3, 4, 5, 6} and X7 = for any values of the 
parameters 71 and 72 in the ranges: 



71 < 
5 



2 

'3' 



< 72 < 0. 



(B24) 
(B25) 



We have not yet found the minimal set of restrictions 
on the parameters that allows the KST constraint evo- 
lution system to be symmetric hyperbolic for any values 
of the characteristic speeds, and for any values of the pa- 
rameter X7- We have limited our numerical experiments 
so far to the regions of this parameter space where both 
the principal and the constraint evolution systems are 
known to be symmetric hyperbolic. 
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